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Abstract. Ensembles of random density matrices determined by various probability 
measures are analysed. A simple and efficient algorithm to generate at random density 
matrices distributed according to the Bures measure is proposed. This procedure may 
serve as an initial step in performing Bayesian approach to quantum state estimation 
based on the Bures prior. We study the distribution of purity of random mixed states. 
The moments of the distribution of purity are determined for quantum states generated 
with respect to the Bures measure. This calculation serves as an exemplary application 
of the "deform-and-study" approach based on ideas of integrability theory. It is shown 
that Painleve equation appeared as a part of the presented theory. 
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1. Introduction 

Random density matrices are a subject of a large current interest. In some cases one 
considers ensembles of random pure states defined on a finite dimensional Hilbert space 
TIk- a natural ensemble is defined by the Fubini-Study measure ^ips^ which is induced 
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by the Haar measure on the unitary group U (K) and invariant with respect to unitary 
rotations. 

In some cases one needs to consider ensembles of mixed quantum states. If the 
dimensionaUty K is a composite number K = MN then an ensemble of random 
mixed states can be obtained by partial trace over an M-dimensional subsystem, 
p = TrM|^)(^|- If random pure states |'^) are distributed according to p^psi one obtains 
in the set of density matrices of order N a family of induced measures [1-3], denoted 
here by iin,m In the symmetric case, M — N, the induced measure is equal to the 
Hilbert-Schmidt measure, which covers the entire set Q, of the density matrices and is 
determined by the HS metric. 

This observation leads to a simple algorithm to generate a Hilbert-Schmidt random 
matrix [2]: a) Take a square complex random matrix A of size N pertaining to the 
Ginibre ensemble [4,5] (with real and imaginary parts of each element being independent 
normal random variables); b) Write down the random matrix 

= IvZ4t' 

which is by construction Hermitian, positive definite and normalised, so it forms a 
legitimate density matrix. Observe that the Ginibre matrix A can be used to represent 
a random pure state of a bipartite system in a product basis, lip) = ^jj^ijlO <8) \ j). 
The above procedure is thus equivalent to setting 

Pus ^ T^n\^){iI^\, where \^l^)&nN®nN (2) 

is a normalized random state taken from the composite Hilbert space of size A^^ according 
to the Fubini-Study measure, while Trjv denotes the partial trace over the second A^- 
dimensional subsystem. 

Another distinguished measure in the space O of quantum mixed states, is induced 
by the Bures metric [6,7], 

Db{p, a) = ^2-2Ti{^pa^pYl\ (3) 

This metric induces the Bures probability distribution, defined by the conditions that 
any ball with respect to the Bures distance of a fixed radius in the space of quantum 
states has the same measure. The Bures metric, related to quantum distinguishability 
[8], plays a key role in analyzing the space of quantum states [9]. The Bures metric is 
known to be the minimal monotone metric [10] and applied to any two diagonal matrices 
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it gives their statistical distance. These unique features of the Bures distance support 
the claim that without any prior knowledge on a certain density matrix acting on space 
Hn, the optimal way to mimic it is to generate the state at random with respect to the 
Bures measure. 

More formally, trying to reconstruct the quantum state out of the results of the 
measurement [11-13], [14, chapt. 3] one can follow the Bayesian mean estimation [15,16]. 
In this approach one starts selecting a prior probability distribution Pq over the set fl of 
all quantum states. Acquiring experimental data one uses them to generate likelihood 
function, multiplies it by the prior and normalizes the result to obtain a posterior 
probability distribution Pi. This distribution reflects the knowledge of an estimator, so 
the best estimation of the quantum state is given by the mean state with respect to this 
distribution, pi = J^pPi{p)dp. If more experimental data arc gathered one continues 
with this procedure to obtain further probability distributions and a sequence of 

expected states, p„ = pPn{p)dp, with n — 2,3, This iterative procedure should 

yield an accurate estimate of the unknown state [17]. 

As the starting point for such a reconstruction procedure one should chose as 
uninformative ("uniform") distribution Pq as possible, so the Bures prior is often used 
for this purpose [18-20]. In practice Bayesian method relays on computing integrals over 
the set fl of quantum states. Since analytical integration is rarely possible, one needs to 
apply some variants of the numerical Monte Carlo method. For this purpose an efficient 
algorithm of generating random states according to a given distribution is required. 
Although the Bures measure was investigated in several recent papers [3,21-25], no 
simple method to generate states with respect to this measure was known. 

Main aim of this work is to solve a few open problems related to the Bures measure. 
We construct the following algorithm to generate random Bures states: a) Take a 
complex random matrix A of size N pertaining to the Ginibre ensemble and a random 
unitary matrix U distributed according to the Haar measure on U{N) [26,27]. b) Write 
down the random matrix 



which is proved to represent a normalized quantum state distributed according to the 
Bures measure. In analogy to the Hilbert-Schmidt case we may also write 



(l + f/)AAt(l + [/t) 



(4) 



Tr[(l + [/)^^t(l + [/t)] 



Pb = 



Trjv|0)((/> 



where 




(5) 
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U G U{N) and {ip) is a random state of a bipartite system used in eq. ([2])- A similar 
construction is also provided to obtain random real Bures matrices. 

The degree of mixture of any state p of size can easily be characterised by its 
purity -P(p) = Trp^. This quantity varies from for the maximally mixed state, 
to unity, characteristic of an arbitrary pure state. Characterisation of purity of 
random states, related to the entanglement of initially pure states before the reduction, 
is a subject of a considerable current interest [29-31]. The average purity is known for 
random states distributed with respect to induced measures, [2,28], and for the Bures 
measure [24] but the distribution of purity is known only for the HS measure for low 
dimensions [31]. For the induced measures the moments of purity were obtained in a 
recent work of Giraud [29]. These results can be rederived by a method involving the 
methods of theory of integrable systems (see [33] and also explanations in part E] of this 
article), which allows to obtain a recurrence relation between moments by deriving 
a differential equation for the corresponding generation function. This differential 
equation is the IV-th Painleve transcendent [35]. Since these moments are already 
known in the literature we will concentrate on a more involved case and derive the 
moments of the purity with respect to Bures measure. Our calculations demonstrate 
practical usefulness of this analytic technique and suggest, it might also be used in 
solving related problems. 

This paper is organised as follows. In section [2] we review probability measures in 
the space of mixed quantum states and provide necessary definitions. In section [3] we 
derive a useful representation of the Bures measure which allows us to construct the 
algorithm based on eq.(j4]). A similar reasoning is provided in section H] for real density 
matrices. In section [5] we analyse the moments of purity for a general class of probability 
measures. Results obtained there are used in section [6] to derive explicit results on the 
moments of purity for Bures random states. Some auxiliary calculations are relegated 
to the appendix. 

2. Ensembles of random density matrices 

We are going to analyse ensembles of random states, for which the probability measure 
has a product form and may be factorised [2,3], 




(6) 
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so the distribution of eigenvalues and eigenvectors are independent. It is natural to 
assume that the eigenvectors are distributed according to the unique, unitarily invariant, 
Haar measure dixu on U (N). Taking this assumption as granted the measure in the space 
of density matrices will be determined by the first factor dz/x describing the distribution 
of eigcnvahies -P(A). 

Consider a class of induced measures iin,m in the space of density matrices of 
size N. To generate a mixed state according to such a measure one may take a random 
bipartite N xM pure state (e.g. an eigenstate of a random Hamiltonian) , and trace 
out the M-dimensional environment. This procedure yields the following probability 
distribution 

dn{p) oc e(p) (5(Trp- l)detp^-^. (7) 

It reflects the properties of density matrices p > and Trp = 1. In the special case 
M = N the term with the determinant is equal to unity and the induced measure 
reduces to the Hilbcrt- Schmidt measure. The matrix p is Hermitian and integrating 
out the eigenvectors of p one reduces dp to the measure on the simplex of eigenvalues 
{Ai, ...Aat} of the density matrix [28], 



dpM{Xi, ...Xn) = Cn,m 1 j n 6>(A^)Af -^dA, 



(8) 



where the squared Vandermonde determinant 



l...Af 
i<j 

appears as a geometric consequence of diagonalisation. The normalisation constant 

Cn M = — Kr~i ^ ^ (10) 

n;=o r(M-j)r(Ar-j + i) 

has been calculated in [2] . 

Furthermore, we analyse the measure induced by the Bures distance, which is 
characterised by the following probability of eigenvalues [3] 

dpBiXr, ..M) = C-l <5 ( ^ A, - 1 j n e{\)Xi'^'dX, if (11) 
The normalisation constant for this measure 
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was obtained in [3, 21] for small and in [23] in the general case. It is easy to see 
that the case = 2 is somewhat special since the denominator in the last factor is 
equal to unity. Incidentally, the Bures measure coincides in this case with the induced 
measure with an unphysical half-integer dimension of the environment, M = 3/2, but 
this observation may ease some computations [25]. 

3. Generating Bures density matrices 

In this section we show that eq. (jlj) may be used to construct an ensemble of random 
states distributed according to the Bures measure. To this end we will rewrite the 
Bures probability distribution corresponding to the measure (flTil in a more suitable 
form, which involves random unitary matrices. 

As a warm-up we shall first consider the induced measure. Let us start with a 
probability measure defined by an integral over random matrices A with respect to the 
Ginibre measure, exp(— TrAA"''), 

Pm(p)cx JdAe--^^^' S^p-^y (13) 

Here A denotes a rectangular complex matrix of dimension A^ x M, and it is assumed 
that M > N. Let us introduce another 5-function by an integral with respect to an 
auxiliary variable s, 

Pm{p) ds J dAe-""'^"^' 5 (^p - 5{s - liAA^). (14) 

After rescaling the matrix variable, A \/sA, the above equation takes the form 
Pm(p) oc P ^e-^s*^^ j dA5{p~ AA^) 5 (l - TiAA^) 

oc e(p) (detp)^^-^ 5(1 - Trp) (15) 

This form is equivalent to (JTj), what proves that random matrices distributed according 
to the induced measure can be generated from rectangular complex matrices of the 
Ginibre ensemble. Taking in particular square N x N matrices one generates random 
Hilbert-Schmidt states according to ([T]). 

To repeat this reasoning for the Bures matrices we will start with a similar ensemble 
defined by a double integral 

PM ocJdAj dHe-^y^^'-"'^^'^ 5 [p - . (16) 
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Here A can be interpreted as a x M Ginibre random matrix, while if is a Hermitian 
matrix of order A^. As in the earher case we introduce a 5-function by integrating over 
an auxihary variable s 

Pb{p) o, ds j dAj rfife-T'-[^^'+^'^-4^] 6 (p - ^) 5(s-TrAAt).(17) 
Rescaling A y/sA leads to 

Pb{p) oc ^e-^s*'^ J dHe-'^'^'p J dA6{p - AA^)6{1 - Tip). (18) 
Performing another rescaling, H H/ y/s, we arrive at 

PM oc r ^e-s'^^^-^'/' U4vll )rXT®(^)^(l - Trp)(detp)^^-^. (19) 



Note that the integration over s gives a constant factor only, which will be absorbed into 
the proportionality relation, while integration over eigenvectors of p gives the squared 
Vandermonde determinant. Furthermore, in the case M = N the last factor equals to 
unity, so the above expression reduces to the Bures measure (fTTll . 

Let us then return to the starting integral f|T6|) and apply another rescaling, 
A -^71=5^. It leads to the following expression 

At this point it is convenient to introduce an unitary variable matrix 

1-iH 



U 



1 + iH 



As shown in lemma 1 proved in Appendix A the 'Cauchy-like' measure dH /[det{l + 
/f^)]^ is equivalent to the Haar measure djiiU) on U{N). Moreover, since 

1 _1 + U 
1 + iH ~ 2 ' 

the above expression is equivalent to 

r r . ( 1+^7/1 4t i+c/t \ 

. MU) I ^Ae-^- S - . (21) 

The factors 1/2 cancel out, so taking a square complex random Ginibre matrix A and 
a random unitary matrix U of the same size we can generate random Bures matrices 
according to the constructive recipe (jl]). Writing a random state {ip) in a product basis, 
\ip) = Xlij^'iK) ® li)) infer that this method of generating random Bures states 
may alternatively be written by eq. ([5]). 
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Figure 1. Distribution P(a) of an eigenvalue a of random matrices of size N = 2 
distributed according to the HS measure (+) and the Bures measure (*). Solid lines 
refer to analytical predictions (|22|) . 



To evaluate advantages of this algorithm in action we have generated in this way 
several random Bures matrices of different sizes. In the one-qubit case, = 2 we 
calculated the distribution of an eigenvalue a = Ai = 1 — A2 for the Hilbert-Schmidt 
and the Bures measures and compared in Fig. [1] our numerical data with analytical 
results obtained in [2], 

PhsW = 12(a - l/2)^ P^(a) = ^^^=1/2!, (22) 

vry a(l — a) 

where a G [0,1]. For larger dimensions we computed the mean traces Km = (Trp™')^ 
averaged over the Bures measure ([H]) . Fig. [2] shows the comparison of the numerical 
data with analytical results following from [3,24], 

, 5N' + 1 3 8N' + 7 

^^'P - 2N{N^ + 2) ' ^^'P - (iV2 + 2)(iV2 + 4)- ^^^^ 

Observe that the average purity of the Bures states is higher than the averages computed 
with respect to the HS measure. This shows that the Bures measure is more concentrated 
in the vicinity of the pure states than the fiat measure. 

Demonstrating practical usefulness of the algorithm to generate random matrices 
according to formula (jll) we may generalise it to get a one-parameter family of 
interpolating ensembles of random matrices. Taking any fixed parameter x from the 
interval [0, 1] and setting y = 1 — x we may construct a random density matrix from a 
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Figure 2. Average traces = {Trp"^)B for a) m = 2 and b) m = 3 for an ensemble 
of Bures random density matrices (x ) and HS random states (+) of size 2, . . . , 10. 
Solid lines represent interpolations of analytical results ([25]) . 



random Ginibre matrix A and a random unitary matrix U, 

^ {yl + xU)AA^iyl + xW) 

Ti[iyl + xU)AA^yl + xW)y ^ ' 

It is clear that for a; = this expression reduces to ([1]) and produces a density matrix 
distributed according to the Hilbert-Schmidt measure, while for x = 1/2 one gets a 
Bures density matrix. Since the Ginibre ensemble is invariant with respect to unitary 
rotations, A UAW, increasing the value of x above 1/2 one interpolates back to 
the HS measure, which is obtained again for a; = 1. Note that the critical parameter 
Xc, at which the transition between both ensembles effectively takes place, is dimension 
dependent, Xc = Xc{N). 

4. Real Bures density matrices 

A similar construction can also be used to construct random real density matrices. To 
generate these matrices according to induced measure [32] 

P^,(p) (X 0(p) 6(1 - Trp)(detp)(*^^^^-i)/2 (25) 

one uses the same formula ([1]) with a random matrix A of the real Ginibre ensemble. 
Note that in the case of real density matrices the Hilbert-Schmidt measure is obtained 
for M = + 1, since in this case the last factor is equal to unity. For this end one needs 
to generate a rectangular real Ginibre matrix A of dimension x (A^ + 1)- 
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To obtain real Bures matrices we begin with an analog of eq. fll6l) in which A is a 
real Ginibre matrix of size x M, while H stands for a real symmetric matrix of size 
iV, 



(26) 



As in the complex case we introduce the 5-function by integrating over s and rescale 
both matrices A and H to obtain expressions 

P|(p) oc ^e-^/"^/2 J ^^e-TtH^P J dA6{p - TiAA^) 6{1 - Tip) 

oc n ^ n XTV®^''^'^^^ ~ Trp)(detp) 



{M-N-l)/2 



This expression coincides with the real Bures measure for M = + 1. In this case we 



perform now another rescaling, A --^==^A. and apply lemma 2 from Appendix A 



In this way we replace an integral over symmetric matrices dH by an integral over 
the measure dpo on symmetric unitary matrices, characteristic of circular orthogonal 
ensemble (COE). The final expression 

\1+u\aa]\1+£!\ 

Pb{p) oc / rfp.(f/) / dAe~^^^' 6 p- „,L., 1 • (27) 



uiN) ■ 'J V TrM^AtMi^ 

allows us to write down the final expression for a real random Bures matrix 

o ^ \1 + V\AA^\1 + V^ 

Tt\1 + V\AAT\1 + V^ ^ 



Here |X| denotes vXJO, while A represents a real rectangular random Ginibre matrix 
of dimension A^ x (A^ + 1), and is a unitary matrix from the ensemble of symmetric 
unitary matrices (COE). To generate such a symmetric matrix one may take any matrix 
U distributed according to the Haar measure on U{N) and set V = UU^ [5]. Also in the 
real case one may design a one parameter ensemble analogous to fl24|) , which interpolates 
between the HS and Bures measures. 

After generating numerically several real random Bures density matrices we 
analysed their spectra. In Fig. [3] we compare the distribution of an eigenvalue P{a) of 
a real one-qubit random state for the HS and Bures measures with the corresponding 
analytical results [23,32], 

Pi«(a) =4|a-l/2|, P|(a) = ^^i£L. (29) 




5. Moment Generating function 

We are going to analyse the moments of a homogeneous function Fq{X) of the eigenvalues 
Aj of degree q for random matrices distributed with respect to the induced measure ([S]) 
and the Bures measure (ITT!) . It is convenient to consider the corresponding Laguerre 
ensembles: 

1...N 

rf/xi,(Ai, ...A^) oc exp(- A.) H - U 0(A.)Af -^dA, (30) 



i<j 



and 



rf/i|(Ai, ...A^) oc exp(- n n 0(A.)Ar'^'dA, . (31) 

. . ~r ^7 

The reason is that the moments and the averages are closely related: 

{F,{\))m = I d/iMF,(A) = f^^^(^.(A))^ (32) 

and likewise 

{F,{X))b = J df,sF,{X) = f^^l^^{F.W)B (33) 

where T{x) is Euler's Gamma function, which can simply be proven by going to spherical 
coordinates. Thus we have the relations for the moments of purity 

(^)M-rp^^^^^(P)M (34) 
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In the same way we may introduce the matrix Laguerre ensembles d/x^(p) and d/x^(p). 
Then we consider the matrix Laplace transforms of these ensembles 

N 

and 



N ^ 

J exp(- Tr Kp)dpi^M = H (1 + K,)^ ^^^^ 



1...N 



which have been calculated elsewhere [24], K is Hermitian K >0 with eigenvalues Ki. 
Prom these we can derive the generating functions for the moments of purity 

N 

^m(^) - J e-^^''^d/.^(p) = e-^W^^)^ 11(1:^^1^.0 ^^^^ 

and similarly 



Applying the matrix differential operator Tr{S/SKy on some invariant function 
it can be expressed in eigenvalues using the Vandermonde determinant A(iif) = 



) (40) 

It is easily seen that this operator is Hermitian. Calculation of all the derivatives, which 
are needed, is not so simple. Instead we make a Hubbard-Stratonovich transformation 
of the exponential operator acting on some invariant function F{K) of degree q, 

e-^-^mK? p(^K) = [ DFexp(-TrF2) exp(2?x/^ Tr F-^)F(li:) 

J 6K K=o 

= j Dyexp(-Try2)F(2iv^y) . (41) 

Here F is a Hermitian matrix. Thus we have reduced this expression to an average 
over the Gaussian unitary ensemble (GUE). We choose the normahsation condition 
Joy exp(— Try^) = 1. For the generating functions we obtain 

N 

= j DY exp(- TY Y^) J] (TT^B^ '''' 
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and similarly 

1...N 

Z^{x) = / Drexp(-Trr2) TT . (43) 

Thus Z^,j{x) is related to some negative moment of the characteristic polynomial in GUE, 
while Zj^{x) is something more complicated - nevertheless also written as some GUE 
average. The above expressions may be considered as a starting point for calculation of 
the moments of the distribution of purity. Since such results were already obtained by 
Giraud for random matrices distributed with respect to HS measure [29], we will not 
discuss this case any further, but we shall rather concentrate on the more complicated 
case of random states distributed with respect to the Bures measure. 

Let us write Zj^{x) as an integral over eigenvalues Yi of Y (with the constant 
Cn = T,N/2Y[^+pY{j) iioi'i^alising the Gaussian measure) 

1...N I...N 



zi,.) ^ c„ / n dK,e- K n (« - y.f n ^^^^ , 
- I n ^=wr ^' n((y^^^ - yrr^))^ . (44, 

Here the Yi integrations run from — oo to +00. Introducing the complex integration 




Figure 4. Contour of integration, C, in eq. ((45|l and its transformation into the contour 
C, right hand side picture, according to the rule Zi 1 + i^/xzi. 



variables Zj = a/1 + 2%^fxYj with Zj^Zj = iy/xdYj we obtain 




where all the Zi run over a contour C starting from e ■ 00 to e"'"*'^/^ ■ 00 passing through 
the saddle Zi = 1 (see figS]). We want to expand Zj^{x) in powers of —x to obtain the 
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averaged moments {P^')^/r\. This means that we have to do saddle-point integration 
for 1/x —>■ oo. The relevant saddle point is = 1. Now we want to expand around the 
saddle point and make the transformation Zi ^ 1 + i^JxZi. 

Z^{x) = Cn [ TT exp(-z2 _ t^zf + xzf/A) ■ A%iz) . (46) 

The new contour C (see figJl]) is such that the integral converges. 

Now we may expand in powers of ^/x, which turns out to become a power series in 
X. In each term we may deform the contour back to the real axis and thus obtain, at 
least for the asymptotic expansion for x ^ 0, the same expansion as for 

Z|(-x) =Cm [ TT dz, exp{-z^ + ^z\ - x4/4) ■ A\{z) . (47) 
with X > 0. Thus 

oo 

Z%{-x) =Y^x^ {P^)%. (48) 
The next section is devoted to the derivation of analytic expressions for moments {P'^)b- 

6. Derivation of moments of purity for Bures measure 

Coefficients of the Taylor expansion of Zj^{—x) in the vicinity of a; = are nothing but, 
up to a constant, moments of traces Trz'^, Trz^ and their powers averaged with respect 
to the probability measure corresponding to GUE. Therefore, our main interest here is 
in calculation of the quantities that below are referred to as Tk^m-, 

T,,^ = ((TrzY(Trz3)2™\ , To,o = 1. (49) 

\ / GUEmxn 

The connection between {P^)^ ^'s is given by the formula 

^ ^ j^y'— m22{m— r)^ 1 

\P)b—/, 7 TTTT^ Tj Tr-mm- (50) 

m=0 ^ J \ J 

See explanations under the formula ( l60i) 

Calculation of ^ for general k and m is a rather nontrivial problem. To do this 
we are going to derive a system of recurrency relations that will allow us, in principle, to 
obtain a closed form of Tk^m for all particular values of k and m. In the basis of derivation 
of such a recurrency lays a so-called "deform-and-study" approach, a string theory 
technique of revealing hidden symmetries. For the first time, this technique as a closed 
calculation method to the problems of Random Matrix Theory appeared in the work [33] 
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by Adler and van Moerbeke where it was utilised to study gap-formation-probability 
integrals over various Unitary Ensembles. Later, this approach was modified to calculate 
distribution properties (such as cumulants) in other random matrix models [34]. 

The celebrated result of the theory of integrable system states that internal 
symmetry of matrix integrals of /? = 2 Dyson's class (encoded in squared Vandermonde 
determinant) leads to highly non-trivial nonlinear relations between combinations of 
averaged traces. One of them which is of primary importance for the approach 
considered is the Kadomtsev-Petviashvili (KP) equation (see [33]). 

The most economic way to work with these relations is to introduce a t-deformation 
into the integration measure. Instead of an original matrix integral 

\...N 

In{x)= Al{z)Y\ dfi{zi;x) (51) 

with some given measure dfi {z; x) depending on a "physical" parameter x, one considers 
the integral depending on an infinite set of auxiliary parameters tk, k = 0,1, . . . 



rivW = ^^^A^(;2)exp|f]4Trz'^|nf^/i(-2.;a;) . (52) 



The KP-equation being written down in the variables t has the following form 



'\+^^-^l^]^ogr^ + Q(^^ogr^ =0. (53) 



In the theory of integrable systems such objects as defined in (1521) are referred to 
as r-functions. Note that tn depends on an infinite set of parameters and satisfies an 
infinite set of relations, one of those is given by ( !53l) . Strictly speaking, one has to 
restrict the number of parameters t^ and then attach correct signs to each of them in 
order to attain convergence of the integral. Since eventually all of them are set to zero, 
we do not need to pay them any special attention. 

The differential equation for the function In{x) is obtained by projecting the KP- 
equation onto the hyperplane t = 0. To perform this projection Adler, Shiota and van 
Moerbeke [33] suggested to use Virasoro constraints (VC) as an additional block giving 
a link between the t^-derivatives in flS5]) and the derivatives over x that supplemented 
the deform-and-study approach to a complete tool applicable to calculation of random 
matrix integrals. On the contrary to the KP-equation, the particular form of the VC as 
well as the way of their derivation is essentially infiuenced by the choice of the measure 
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(i/i {z\ x). The basic idea of derivation can be expressed in the following simple form 

eO, g = -l,0,l,... 



5e 



while the details may vary from one model to the other, and they can be demonstrated 
much simpler by studying of our particular examples. 

The final step of the approach is to resolve the obtained VC and KP-equation jointly 
on the hyperplane t = to bring a closed equation for 2^{x). Substitution of T^^x) in 
the form of a Taylor series into the obtained equation gives rise to a recurrence relation 
for the coefficients of the series and subsequently for the sought moments T^. eq. (H9|) . 

Below we show how this approach can be applied to calculate the moments ^ 
defined by ( 149|) . All details of calculations are given in two appendices ( Appendix B 



and [Appendix C| ), below we give only a plan of the calculation program. 

We start with the derivation of the recurrence relation for T^ q- This choice is 
dictated by two reasons, first, by the relative simplicity of this case and, second, 
because the moments Tkfi serve as an initial condition for calculation of the higher 
order moments. 

To derive expressions for the moments T^^q = \(Trz'^)'^\ we consider an 

auxiliary integral, Jn{,x), 

2N{N-l)/2 I- l---^ 

Obviously, the sought moments can be obtained as coefficients in the expansion of this 

integral into the Taylor series in the vicinity of x = 0. Thus, deriving a differential 

equation on Jn{x) by using the method discussed in the main section enables us to link 

the moments by a recurrence relation. 

Direct application of the "deform-and-study" approach to the integral ( l54l) does not 

help in the derivation of a differential equation for the function Jn{x) (this is discussed 

in details in [Appendix Bj ). However, due to its symmetry the integral Jn{,x) can be 

represented product of two simpler integrals: 

2N{N-i)/2 [s+(a;)S-(a;), N = 2k- 

= ,,,,^jv < (55) 



where 




N = 2k + l 



^kix) = TT A^(^) ^''e-^^-^^^dz,, V = ±1. (56) 
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Integral can be investigated with the help of the announced approach. Details of 



calculations and results are given in [Appendix B[ Here, however, we make an important 



remark. There is another way to derive equation ( IB. 51) from the appendix. It is a 
consequence of the known result by P.J. Forrester and N.S. Witte [35]. They showed 
by using other methods that the matrix integral (here we use their original notation for 
integral) 

En{s- y) = -, Aliz) His - z.fe-^Uz, 

^- J(-oo,s]" 

is expressed in terms of a solution of the Painleve IV equation. Namely, 

-f log^„(s;z/) = 0(s), (57) 
as 

where 0(s) satisfies the differential equation, which is the Painleve IV equation written 
down in the Chazy form 

(p'" + 6{(p'Y + 4[2(n - i^) - s^](p' + 4s0 - 8nu = 0. (58) 

Then, due to the relation 



' 2^ 

one can restore ( ]B.5p from ( 158|) and ( l57|) by appropriate change of variables. In spite of 



the fundamental character of the obtained Painleve, the equation (1B.5P derived by our 



regular method is more convenient for the analysis that is done in [Appendix B 

To derive the recurrence relation for the Tk^m with the help of "deform-and-study" 
approach as well as in the previous case we define an auxiliary integral 

nN{N-l)/2 r ^ 

^) = ^/.rr^+ip,, / ^Niz) n e-^^-^-.^rf., (59) 
^ ^ Uj=i r(j) Jn^ fJl 

To make the further procedure of derivation successful we have introduced the extra 

parameter y. Appearance of more then one variable apparently lead to a differential 

equation in partial derivatives in both variables. 

The form of integral ( l59l) implies that one can seek the solution of the obtained 

equation in the form of a series in x and y: 

CO m k 

J^ix,y)= (60) 

k,m=0 

Expansion (!60!) is used to derive relation (l50l) between averaged moments {P^)b and 
Tkrn- Indeed, it is enough to compare coefficients standing at equal powers of x in the 
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Taylor expansion of the both sides of the identity J]^{x/A, y/x) = Z^{—x). Expansion 
in the left hand side can be obtained from eq. fl60l) after the change x ^ x/A; y —>■ ^/x, 
while in the right hand side one uses the Taylor expansion ( l48l) . 

Explicit expressions for the first several moments ,„ are given in 



Appendix C 



The higher moments of purity follow from the general expression (!35|) which being 
combined with the relation (1501) gives 

T{N^/2) ^ (_i)'-'»22("^-'-)r!, 



((Trp^)Oi. 



r(A^V2 + 2r)^ (r-m)!(2m)! 



-T, 



r—m m- 



(61) 



To demonstrate the ability of our approach we calculated the first three moments fir 
explicitly. The first moment fii coincides with the mean trace (Trp^)^ given by eq. (123|) . 
Expressions for other two moments are reproduced below 

5(5iV^ + 47Ar2 + 32) 



1^2 



((Trp^)')B 



/is = ((Trp2)') 



B 



5{25N^ + 690N^ + 6015A^^ + 8750 A^^ ^ 1^52) 



(62) 



8A^(A^2 + 2)(A^2 ^ 4) ■ ■ ■ {N^ + 10) 
As shown in Fig. [5] our numerical data obtained by generating random Bures states 
according to the method presented in Sec. [3] coincide with the above analytical results. 




Figure 5. Average moments of purity /ifc = ((Trp'^)''')^ with k — 1,2,3 for an 
ensemble of Bures random density matrices of size N — 2, . . . , 10. Solid lines represent 
interpolations of analytical results (|62p . 
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In this work we proposed an explicit construction to generate random density matrices 
according to the Bures measure. A single complex random Bures density matrix of size 
N is obtained directly out of one complex random square matrix A from the Ginibre 
ensemble and one random Haar unitary matrix U of size N. Similarly, to generate a real 
Bures state of size A^, it is sufficient to have an rectangular, N x (A^ + 1) real Ginibre 
matrix and a random symmetric unitary matrix V. These practical recipes are not only 
simple but also economic and allow one to form the Bures prior, useful as an initial step 
to apply the quantum Bayesian approach. 

Studying the distribution of random Bures states we have analytically determined 
the moments of their purity. These results, derived by means of the theory of integrable 
systems, reveal the power and usefulness of this analytic technique. Note, that the 
Painlcvc IV transcendent appears in a natural way as an intermediate step of our 
calculation. 
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Appendix A. Measures on Hermitian and uniteiry matrices 

In this section we prove two lemmas which allow us to replace integration over the set 
of Hermitian matrices by integration over unitary matrices with respect to the circular 
unitary ensemble (CUE) and circular orthogonal ensemble (COE), respectively. 

Lemma 1. The measure dH/[det{l + H'^)]^ on the set of Hermitian matrices 
H — of order N, is equivalent to the Haar measure on U{N), which corresponds to 
CUE. 
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Proof. We start with a substitution, U = {1 - iH)/{l + iH) = 2/(1 + iH) - 2. 
Thus the measure dU = — ^ ^ rr idH ^ imphes that U^^dU = '^.„ idH ^ This 
allows us to write down an explicit expression for the metric 

ids)' = TT{U^'dU){U-'dUy = ATij^^dHj^dH, (A.l) 

which implies the measure 

di^(U) = (A.2) 

^ [det(i^)]^ ^ ^ 

and completes the proof. 

A similar lemma can be formulated for the measure on U{N) related to COE. 

Lemma 2. The measure dH/[det{l + if 2^j{Af+i)/2 ^.j^^ gyjjijjietric matrices 
H = = W of order N, is equivalent to the COE measure fio on U{N). 

Its proof is analogous to the previous one and the exponent (A^ + l)/2 is related to 
the number of A^(A^+ l)/2 of the independent variables of a random symmetric matrix. 

Appendix B. Recurrence relations for T^^ 

Proof of ( 1551) . To prove fl55l) we notice that the matrix integral Jn{x) has a 
determinant al representation: 



Jn{x) ~ det 

0<j,i<n-l 



Due to the symmetry of the function under integral the matrix of moments in the above 
determinant has a chessboard structure with zeros on all "white squares", i.e. for the 
elements with i + j = 2k + 1 {k = 0,1, . . .). This type of determinants can be reduced 
by permutation of lines and rows to a determinant of block-diagonal matrix and as a 
result to a product of two determinants. In our particular case it gives 

n^-%i 1 [^*+^'] = , [^2i+2i+2] det [^l2^+2j] = det ■] dct .] 

0<''^<^-^ 0<i,i<[fj 0<i,i<[f] 0<i,j<[f\ 0<i,i<[f] 

where [•] and [•] denote the integer part of a real number, so that [a;J < a; < [a;] 
and the momentum matrix fif^j is obtained from H2i+2j and ^2i+2j+2 by the change of 
variables z'^ ^ z in corresponding integrals, 

jjf^^{x) = I z'+^^'/^e-'-''''dz 
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Now returning back from the determinants to the matrix integral representation we 
obtain fl55l). 



r-function and Virasoro constraints. The integral corresponds to the 

r-function of the form 



f,{t} = l / A,(^)n^f exp 

■ "^^+ .7 = 1 



-Zj - x^l + ^ t^zj 
1=1 



dzj , (B.l) 



here, the parameter u stands for ±1. 

To derive the Virasoro constraints (VC) first one has to choose an appropriate 
change of variables. The general recipe says that this transformation must be chosen as 



Zj t-^ Zj 



+ 5ezf-^f{z,) \{{z, - afc) , g = -1, 0, 1 . . . 

k 

where are the boundaries of integration domain without both infinities (if they are 
presented) and with excluded zeros of the polynomial function f{z). The function f{z) 
is, in turn, related to the original integration measure through the parametrisation 



with h{z) = ^ Qkz'' , f{z) = J2f> 



dz fiz) 

^ ' fc=0 fc=0 

where Vi^z^ is a confinement potential, in our case 

dV{z) 2z + Axz'^-iy 



kZ 



V{z) = z + xz — — log z 



dz 2z 



and correspondingly f{z) = 2z. Since the only zero of f{z) coincides with the only finite 
boundary point of the integration domain, z = 0, we should use the shift of the form 

Zj Zj + 25ez]^^ g = -1, 0, 1 . . . (B.2) 

Substitution of this change of variables into the integral fIB.ip and variation over 6e give 
rise to an infinite number of VC 

2 T^n^T + 2 ^''^rnjr. + V— 2— Ax— = 

m=0 '^'''"'^ m=l ^'^q+ra C'^g+l OLq^2 C^g+S 

where the first term originated from the squared Vandermonde determinant and the 
volume element 11^=1 '^^ji second term corresponds to i- deformation, and the other 
three are the contributions from the measure. Note that the operation of differentiation 
over to reduces to multiplication by the integral dimension, ^ = k. 
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Projection of KP onto the hyperplane t = 0. To perform the projection of the 
KP-equation fl53|) onto the hyperplane t = one needs to know the following derivatives 



log ffc 



dt\ 



d"^ log Tk 



dtj 



log Tk 



t=0 



dtidt3 



t=0 



the expression for the second derivative over t2 immediately follows from the observation 

dfk 



dx 



dtn 



The same observation allows us to rewrite the first two VC (g = — 1 and g = 0) in the 
form that involves among t-derivatives also x-derivatives (below we use the notation 
g{x;t) = log ffc {t}) 



2k'^ + 2^ mtr. 



m=2 



-^ + 2t,^ + ku-2^ + Ax^ = 0- 

Olm Oil Oil OX 



2t — + V— + 2— - Ax— 

dx dti dx dt^ 



0. 



(B.3) 
(B.4) 



dti dtm+i 
These two equations give all necessary information. Indeed, from (1B.3I) we obtain 
derivatives over ti: 

dg{x; t) 

t=o 



dti 




d'^g{x] 


t) 


dti 




d^g{x; 


t) 


dtf 



k{k + u) dg{x;0 
z 1- 2x 



k(2k + v) 



+ 2x 



dx 
dg{x- 0) 



t=o 



t=o 



?,k{2k + z/) + 12a; 



dx 
dg{x; 0) 



+ 4(x^) ~g{x;0) 



dx 



+ 44 x 



dx 



+ 48 X 



d_ 

dx 



d 



g{x;0) 

4 



^(x;0) + 16(x— ) gix;0). 



Then from flB.4l) one can get the mixture derivative over ti and 



d^g{x; t) 



dtidts 



t=Q 



1 

4x 



(4fc + z/) 



d^g{x; t) 



dti 



d dg{x;t) 



t=o 



, dg{x; 0) 
dx 



Substitution of these results into the KP-equation gives rise to a nonlinear equation in 
partial derivatives of the function g{x;0): 

1 3 
- -k {8k'^ + Qku + + -k{2k + u) [2k'^ + kv + 2)x- 

2(1 + 3(4A; + z/)x - 6 {Qk'^ + ?,kv + lO) x^) g' 

-X (1 + 4(4A; + u)x - 12(4P + 2ku + 2^)x^) g" 

+216x^ {g f + 288x^^'^" + 96x^ {g" f + 144x^^(^) + IGx^^^^^ = (B.5) 
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Note that the procedure of joint resolving of KP-equation and VC's fails if we try 
to apply "deform-and-study" approach directly to the integral fl54l) . In this case the 
g-th VC contains the term -Qf^, which becomes already at g = —1 a derivative over ^3. 
This gap in derivatives makes the system KP-VC unresolvable. 

Recurrence relation and some explicit results for Tfco- To derive the 
recurrence relation for the moments we substitute, first, the Taylor expansion of 
the function g'^{x) (equipped with an extra index z/ = ±1): 

g^{x) = logS^(x) , El{x) = al V(-l)'"c^™ — , 

ml 

m=0 

where 

k-l 



(B.6) 



3=0 



27' 



into equation (IB.Sp . In the derivation we used a general relation 



X 



{-xy 



0=0 



dx^ 



A A ■ 

?! ^ — ' m! ? — m — c ! 

j=e m=o ' 

As the result we obtain the recurrence relation for the coefficients (it is assumed that 
the summation up to a negative limit is an identical zero) c^^ (m > 1) 



(2 + m)~cl ,~cl = -k{2k + v) {Ak + z/) ^ 



3=0 



k j^k m—j 



m—1 



+^E 

j=0 



m—1 
j 



-k{2k + u){2e + ku + 2)dl^cl^_^_, + 6{4k + u)~clfcl 



J + 1 



1 L^I- 



k j+l^k m—j 



m-2 

+4j(m- 1)^ 

j=0 



m — 2 
j 



m—3 



+ ^ku + n~circl - [Ak + z/)cl 



+{Ak + vYcirc 



j+l k m—j — 1 



j^km-j ~ _|_ ^-^^k j+l^k m-j 



'12j{m- l)(m-2)^ 

j=0 



m — 3 
J 

m—A 



{Ae + 2ku + 7)clj^,cl^_j_, 



-{AP + 2ku + 25)cljS[ 



k m—j — 1 



-lUj{m - l)(m - 2)(m - 3) ^ 

j=0 



m — A 
j 



^k j+l^k m—j— 2 ^k j^k m—j—1 
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m—5 
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+16j(m - l)(m - 2)(m - 3)(m - 4) ^ 

j=0 



m — 5 
j 



'^^k j+2^k m-j-3 ^^k j+l^k m-j-2 



(B.7) 



■"k j k m—j — 1 ■ 

It can be resolved with the only initial condition c^q = 1 that follows from the particular 
choice of the parameter a'^ in the expansion flB.6l) and the definition fl56|) of 

The sought moments o one can find from (1551) by comparison of the the coefficients 
of Taylor expansions in both sides: 



£=0 



2Ar(Ar-l)/2 

vr^/^nf=ir(j) 



£=0 m=0 



^! ^ \m 

e=o m=0 ^ 

Then, after a slight massage we arrive to a simple formula 



^k m^k £-mi 



^k m^k+1 l-mi 



N = 2k: 



N = 2k + l 



T,o = ((Tr^Y) =E 

\ / GUEjSfxN r. 

m=u 

Below we reproduce the first three moments 



^k m^s l—m ' 



k=[N\ , s= \N] 



(B.8) 



^0,0 — 1; 

_ AT 
Ti,o - ^ + 

Tsfi = —r- + —r. — + + -^TT- + —■ 



64 



32 



16 



Appendix C. Recurrence relations for T; 



k,m 



r-function and Virasoro constraints. To define the r-function for this case we 
perform t-deformation of the measure in the original integral (159|1 . so that 



TNit} 



1...N 



A%{z) H exp 



k=l 



dzj, (C.l) 



The extra dependence of r on one more additional parameter allows to resolve 
successfully the KP-equation and the VC. A similar question was discussed in details in 
the paragraph under equation flB.51) . 
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To derive the VC for the r-function (IC.ip we use the transformation of the form 



g = -l,0,l 



then the general form of VC reads as follows 



dtmdtm^q ^ <9tg+m dtq^2 

m=0 ^ m=l ^ ^ 



dt 



q+3 



dt 



g+4 



Observing that 



dy dt3 
dx dt4 ' 

we can rewrite the first two VC (g = —1 and g = 0) in the form {g{x, y; t) = logrAr {t}) 

dg dg ^ dg dg 



Nt, + V mtm^r^ - 2-f + 3y-f - 4a;- 
^ dtm-i dtl ot2 oy 



0; 



m=2 

oo 



/\r2 , \^ + , ^ dg dg 

N + } mt^— 2— + 3y— + Ax— = 0. 

^ dtm ot2 oy dx 



(C.2) 
(C.3) 



Projection of KP onto the hyperplane t = O.To perform the projection of the 
KP-equation (1531) onto the hyperplane t one needs to know the following derivatives 



d'^g{x,y;t) 



dt\ 



d'^g{x,y;t) 



t=o 



dtidy 



d^g{x,y]t) 



dt^. 



and 



dtf 



t=0 



t=0 ^^2 

The derivative over t2 can be expressed from fIC.SP : substitution of the latter into f lC.2p 
helps to find the derivative then the necessary projections can be found by 
subsequent differentiations over ti and t2. As the result we obtain 



d'^g{x,y; 


t) 




dtidy 




t=o 


d^g{x,y; 


t) 




dtl 




t=o 


d^g{x,y; 


t) 




dtl 




t=0 


d^g{x,y; 


t) 





I (63^^ - 88x) ^(°'^) + ^ (8x - + ^ - 8x) 



,(1,1). 



9tf 



t=o 



16 
45 
256 



(5696x^1/ - 14544x1/^ + 6237y^) ^ 



5^ n(o,i) 
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/10287xV _ 57429x1/'^ 331695i/6 _ 2) 
V 4 16 256 ^ )^ 

+—xy {U2Ax'^ - 11988x?/2 + 6939?/^) ^^^'^^ 
16 

-^x (8x - 9?/2) (32x2 - 396xy2 + 297?/^) ^^^'^^ 

+ 1^^ (8^ - V)' (631/2 _ 88^) ^(0,3) _ (8x - 9y2)^^(i-3) 

+ ^x2y2 (8^ _ ^y2Yg%2) + ^ (8x - 9y2)%(o.^). 

Above it is assumed that in the right hand side the function g is taken at t = 0. The 
standard notation for partial derivatives, g'-'^'"^' = gk^ gmy ; is also used. 

Substituting these terms into KP, eq. flS^ . we arrive to a nonlinear equation in 
partial derivatives on the function log J'n{x,y) with the maximal derivative of the forth 
order and quadratic nonlinear terms. This equation can be rewritten in the form of 
equation on the function Jn{x,y) itself. The explicit form of this nonlinear (all terms 
are quadratic in Jn{x,y)) equation is too cumbersome to be reproduced on paper. 

Recurrence relation and some explicit results for T/^ m. Substitution of the 
expansion (HI]) . 

fc,m=0 

gives rise to the sought recurrence relations for the coefficients Tk^m- They are easier to 
be handled by using computer rather then a pencil, and here, we do not reproduce them 
in any form, all necessary calculations were done with the help of computer. Below we 
give results for the first several moments. 



_ 3N 3N^ 

To,i - X + 

471A^2 225A^^ 3A^^ 
TiA = + -—— + 



32 16 4 ' 

_ 4563Ar2 675Ar^ 27N^ 

495iV 82335iV3 8673A^5 555Ar^ 3N^ 

^'^ 4 128 32 32 8 ' 

25515A^ 1194939A^3 292383A^^ 1323A^^ 27A^9 

^'^ 32 256 128 8 8 ' 

382725A^ 19566765A^3 2713095A^5 59535A^^ 405A^9 

64 512 128 32 8 
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